rm(list=ls())
setwd("/Volumes/Disk2/Future_PM/Dataverse")
library(fields); library(maps); library(ncdf4);library(abind)
open.ncdf=nc_open;get.var.ncdf=ncvar_get;close.ncdf=nc_close

source('Function/get_geo.R')
source('Function/get_met.R')
source('Function/functions.R')
source('Function/filled.contour3.R')

#================= 05x0667 ================
filename='Data/SO4_2004-2012.nc'
datafile = open.ncdf(filename)
lon=seq(-180,177.5,by=2.5)
lat=seq(-90,90,by=2)
ind1=which(lon>=-130 & lon<=-65)
ind2=which(lat>=25 & lat<=50)
lon.out=lon[ind1]
lat.out=lat[ind2]

time=1:27
T=get.var.ncdf(datafile, start=c(ind1[1],ind2[1],1),count=c(length(ind1),length(ind2),length(time)),varid='total_T')
PBL =get.var.ncdf(datafile, start=c(ind1[1],ind2[1],1),count=c(length(ind1),length(ind2),length(time)),varid='total_PBL')
cloud =get.var.ncdf(datafile, start=c(ind1[1],ind2[1],1),count=c(length(ind1),length(ind2),length(time)),varid='total_cloud')
PBLcloud =get.var.ncdf(datafile, start=c(ind1[1],ind2[1],1),count=c(length(ind1),length(ind2),length(time)),varid='PBL_cloud')

PSO4_gas =get.var.ncdf(datafile, start=c(ind1[1],ind2[1],1),count=c(length(ind1),length(ind2),length(time)),varid='PSO4_gas')
PSO4_O3 =get.var.ncdf(datafile, start=c(ind1[1],ind2[1],1),count=c(length(ind1),length(ind2),length(time)),varid='PSO4_O3')
PSO4_H2O2 =get.var.ncdf(datafile, start=c(ind1[1],ind2[1],1),count=c(length(ind1),length(ind2),length(time)),varid='PSO4_H2O2')

OH=get.var.ncdf(datafile, start=c(ind1[1],ind2[1],1),count=c(length(ind1),length(ind2),length(time)),varid='total_OH')
O3 =get.var.ncdf(datafile, start=c(ind1[1],ind2[1],1),count=c(length(ind1),length(ind2),length(time)),varid='total_O3')
SO2 =get.var.ncdf(datafile, start=c(ind1[1],ind2[1],1),count=c(length(ind1),length(ind2),length(time)),varid='total_SO2')
H2O2 =get.var.ncdf(datafile, start=c(ind1[1],ind2[1],1),count=c(length(ind1),length(ind2),length(time)),varid='total_H2O2')
close.ncdf(datafile)
T[T<=273]=NA



#===============================
# T  PBL  cloud   PBLcloud  
# PSO4_gas  PSO4_O3  PSO4_H2O2 
# OH O3   SO2   H2O2

date=cbind(rep(2004:2012,each=3),rep(6:8,9))
for(month in 6:8){
	ind=(date[,2]==month)
	for (i in 1:length(lon.out)){
		for (j in 1:length(lat.out)){
			T[i,j,ind]= mov.detrend(T[i,j,ind],len=5)
			PBL[i,j,ind]= mov.detrend(PBL[i,j,ind],len=5)			
			cloud[i,j,ind]= mov.detrend(cloud[i,j,ind],len=5)	
			PBLcloud[i,j,ind]= mov.detrend(PBLcloud[i,j,ind],len=5)	
			PSO4_gas[i,j,ind]= mov.detrend(PSO4_gas[i,j,ind],len=5)	
			PSO4_O3[i,j,ind]= mov.detrend(PSO4_O3[i,j,ind],len=5)												
			PSO4_H2O2[i,j,ind]= mov.detrend(PSO4_H2O2[i,j,ind],len=5)	
			OH[i,j,ind]= mov.detrend(OH[i,j,ind],len=5)												
			O3[i,j,ind]= mov.detrend(O3[i,j,ind],len=5)	
			SO2[i,j,ind]= mov.detrend(SO2[i,j,ind],len=5)												
			H2O2[i,j,ind]= mov.detrend(H2O2[i,j,ind],len=5)												

		}
	}
}


cal.slope=function(x,y){
	slope=array(NA,c(dim(x)[1],dim(x)[2]))
	sig=array(NA,c(dim(x)[1],dim(x)[2]))
	for(i in 1:dim(x)[1]){
		for(j in 1:dim(x)[2]){
			fit<-lm(y[i,j,]~x[i,j,])
			slope[i,j]=coef(fit)[2]
			sig[i,j]=summary(fit)[[4]][2,4]
		}
	}
	slope[sig>0.05]=NA
	return(slope)
}

plot.slope=function(spdata){
	ind=(date[,2]>=6 & date[,2]<=8)
	cr=find.Rlim(0.05,sum(ind))
	slope= cal.slope( month_T[,,ind],spdata[,,ind])
	slope[!land.mask]=NA
	plot.field(slope, lon.out, lat.out,type='sign',mai=c(0.1,0.1,0.15,0.1),legend.mar=4.5,legend.width=1.3)
	clines=contourLines(lon.out, lat.out, slope,level=c(1))
	if(length(clines)>0)for(k in 1:length(clines)){lines(clines[[k]]$x,clines[[k]]$y,lwd=1,col=1,lty=2)}
}


ss=load("Data/land_mask.Rdata")
land.mask=sp.dissolve(mask,lon.in-360,lat.in, lon.out,lat.out)>0.1

#=========================
mai=c(0.1, 0.1, 0.2, 0.1)
mgp=c(1.2, 0.4, 0)
tcl=-0.2
ps=12
close.screen(all.screens = TRUE)

m <- rbind(c(0,1,0.9,1),
	c(0,0.33,0.6,0.9), c(0.33, 0.67, 0.6, 0.9), c(0.67, 1, 0.6, 0.9),
	c(0.15,0.85,0.24+0.3,0.54+0.3)
)
split.screen(m)

screen(1)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
text(x=0.5,y=0.2,"Slopes of sulfate production rates with temperature",cex=1.1)

screen(2)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
ind=(date[,2]>=6 & date[,2]<=8)
slope= cal.slope( T[,,ind], PSO4_gas[,,ind])/10000
slope[slope>=10]=10;slope[slope<=-10]=-10
slope[!land.mask]=NA
image(lon.out,lat.out, slope,zlim=c(-10,10),col=rwb.colors(32),xaxt='n',yaxt='n',xlab='',ylab='');
map("world", add = TRUE)
box()
title('(a) via OH (gas)',cex.main=0.9,font.main=1)

screen(3)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
ind=(date[,2]>=6 & date[,2]<=8)
slope= cal.slope( T[,,ind], PSO4_H2O2[,,ind])/10000
slope[!land.mask]=NA
slope[slope>=10]=10;slope[slope<=-10]=-10
image(lon.out,lat.out, slope,zlim=c(-10,10),col=rwb.colors(32),xaxt='n',yaxt='n',xlab='',ylab='');
map("world", add = TRUE)
box()
#title('(b) via H2O2 (aqueous)',cex.main=0.9,font.main=1)
title(bquote('(b) via'~H[2]*O[2]*' (aqueous)'),cex.main=0.9,font.main=1)

screen(4)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
ind=(date[,2]>=6 & date[,2]<=8)
slope= cal.slope( T[,,ind], PSO4_O3[,,ind])/10000
slope[slope>=10]=10;slope[slope<=-10]=-10
slope[!land.mask]=NA
image(lon.out,lat.out, slope,zlim=c(-10,10),col=rwb.colors(32),xaxt='n',yaxt='n',xlab='',ylab='');
map("world", add = TRUE)
box()
title(bquote('(c) via'~O[3]*' (aqueous)'),cex.main=0.9,font.main=1)

screen(5)
 mai=c(0.1, 0.1, 0.05, 0.2)
 par(mai=mai, mgp=mgp, tcl=tcl, ps=9)
image.plot(slope,legend.shrink=1,legend.only=TRUE,zlim=c(-10,10),col=rwb.colors(32),horizontal=TRUE,legend.width=2.5,axis.args = list(mgp = c(0, 0.5, 0),tcl=-0.2,padj=-1
))
text(bquote(''~10^4~kg~K^-1*''),x=par("usr")[2]-0.15,y=par("usr")[3]+0.625,srt=0, adj = 0,xpd=TRUE,cex=1.3)

